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We present a minimal model for the formation and migration of aeolian sand dunes. It combines 
a perturbative description of the turbulent wind velocity field above the dune with a continuum 
saltation model that allows for saturation transients in the sand flux. The latter are shown to 
provide the characteristic length scale. The model can explain the origin of important features 
of dunes, such as the formation of a slip face, the broken scale invariance, and the existence of a 
minimum dune size. It also predicts the longitudinal shape and aspect ratio of dunes and heaps, their 
migration velocity and shape relaxation dynamics. Although the minimal model employs non-local 
expressions for the wind shear stress as well as for the sand flux, it is simple enough to serve as a 
very efficient tool for analytical and numerical investigations and to open up the way to simulations 
of large scale desert topographies. 

PACS numbers: 45.70.Mg, 45.70.Qj, 47.27.-1, 51.10.+y 



I. INTRODUCTION 



Sand dunes develop wherever sand is exposed to an 
agitating medium such as air or water that lifts grains 
from the ground and entrains them into a surface flow. 
The diverse conditions of wind and of sand supply in 
different regions on Earth give rise to a large variety of 
shapes of aeolian dunes |[ |^. Moreover, dunes have 
been found on the sea-bottom and even on Mars ^, |j. 
Despite the long history of the subject, the underlying 
physical mechanisms of dune formation are still not very 
well understood. How are aerodynamics (hydrodynam- 
ics) and the particular properties of granular matter act- 
ing together to create dunes? How is the shape of a dune 
maintained when it moves? Since the macroscopic phe- 
nomena of interest are separated by many orders of mag- 
nitude from the grain scale and involve various coupled 
nonlinear processes such as turbulent air flow and grain 
hopping ("saltation"), one is bound to devise some sim- 
plified models in order to address such questions. We will 
argue that approximate numerical models can only be 
successful if based on a sound qualitative understanding 
of the problem. Therefore, our main aim is to identify the 
key mechanisms underlying dune formation and migra- 
tion and incorporate them into a working minimal model 
of aeolian sand dunes, and we will emphasize generic as- 
pects over the more specific details. For dcfiniteness, the 
reader may find it helpful to think of isolated transverse 
dunes or crescent-shaped barchan dunes as major appli- 
cations of the model. A schematic sketch of the height 
profile of a barchan is shown in Fig. ||. The broad phe- 
nomenology of aeolian and submarine land forms pro- 
vides a large number of different characteristic structures 
that can certainly not all be described by the same simple 
model developed with the specific examples of barchan or 
transverse dunes in mind. However, we expect that our 
approach is amenable to future adaptations that make it 



applicable to a broader class of sand topographies on the 
one hand, and for quantitative investigations of more spe- 
cific questions on the other hand. Although the minimal 
model refers only to rather generic properties of the wind 
velocity field and the laws of aeolian sand transport, it 
can make interesting predictions about the surface pro- 
file, the development and position of the slip face, dune 
migration etc. that are insensitive to the simplifying as- 
sumptions. The main features of the model were already 
briefly presented in a recent Letter |^ . The present con- 
tribution gives a more comprehensive discussion of the 
model and tries to communicate its precise definition as 
well as its major predictions to an interdisciplinary read- 
ership. The model, as presented here, is restricted to 
a two-dimensional {2d) slice of a dune parallel to the 
unidirectional wind. (A generalization to Zd problems 
is in preparation.) A further restriction is the neglect 
of ripples and direct slope effects onto the sand trans- 
port outside slip faces. Although they have successfully 
been incorporated into continuum sand transport mod- 
els [Q, |[ || similar to our own , we chose to disregard 
them for the present purpose and leave their integration 
to future work. 

The paper is organized as follows. In the next intro- 
ductory section we summarize some background knowl- 
edge and basic definitions. We will also introduce a naive 
"zeroth order" description of the wind shear stress and 
the induced aeolian sediment transport. Its instructive 
failure to produce dune-like steady-state solutions will 
be a guide for identifying two relatively small effects (the 
upwind shift of the maximum of the shear stress with re- 
spect to the topography and the saturation transients in 
the sand flux) as key ingredients of a proper description 
of structure formation by aeolian sand transport. We 
will moreover derive the scaling behavior of the migra- 
tion velocity for translation invariant heaps and dunes 
of different size but similar shape based on very general 
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FIG. 1: Sketch of a barchan dune. Sand is eroded by the wind 
on the upwind or "stoss" side and transported to the brink. 
Strong deposition occurs due to flow separation behind the 
brink. On the downwind or "lee" side, sand slides down at 
the angle of repose (about 32° — 35°) over a concave slip face. 



grounds. Sections [II and |^ are devoted to the defini- 
tion of the minimal model, i.e., to the modeling of the 
air shear stress exerted onto a heap of sand and the in- 
duced sand transport, respectively. The first step builds 
on turbulent boundary layer calculations developed in 
a series of publications mainly by Hunt and coworkers 
PI , [l2| [l^ [l^ , [l^ , the second one on a previous con- 
tribution |10| by the present authors. Only the most 
pertinent results of these earlier developments will be 
summarized here. In the remainder, we will derive some 
important predictions of the model for the central slice 
of a barchan dune or transverse dune. In particular, we 
will demonstrate that there is a minimal dune size. Al- 
though we will thereby gain interesting results, these are 
rather meant to be illustrative examples of possible ap- 
plications of the model. By no means do we attempt 
to provide a complete analysis of its predictions, and it 
should become obvious that much more remains still to 
be done. Finally we will summarize our main results and 
speculate about probable consequences of the present 2d 
theory for 3d topographies. 



II. GENERAL 
A. Aeolian sand transport 

Before going into the description of the model, we want 
to recall some general background and to introduce some 
quantities of major interest. First of all, for convenience, 
we will usually refer to dunes without slip face as heaps. 
Further, we will sometimes find it helpful to focus on 
isolated heaps or dunes on bedrock, although most of 
our discussion is not restricted to this situation. 

The key quantity for the description of the formation 
and migration of sand dunes and heaps is the local hor- 
izontal surface velocity v{x,t) of a sand height profile 
h(x, t) at all positions x and times t. Via mass conserva- 
tion it can be related to the erosion rate ^q{x, t) (nega- 



tive erosion is deposition), where the sand flux q{x,t) is 
defined as the mass of sand transported per unit of time 
across a hyper-plane transverse to the wind direction. 
More precisely, since we want to specialize our discussion 
to a 2d slice parallel to the unidirectional wind veloc- 
ity, the hyperplane is a vertical line and g is a flux per 
unit width. Mass conservation then takes the form of a 
continuity equation for the height profile 



dh{x,t) dq{x,t) 



dt 



dx 



(1) 



with Qg the density of the sand bed. 

With Eq.(Q) one can write the position dependent mi- 
gration velocity at a given time t as 



v[x) = Qs 



(2) 



where we have introduced the shorthand notation 
f'{x) = df{x)/dx. At this stage we can already get some 
physical insight by observing that this equation needs 
special attention at the top of a heap or dune, where we 
expect the denominator to vanish. For v to remain fi- 
nite at the crest as required in the steady state, there 
are in general only two possibilities. Either the sand flux 
q is fine tuned so that the erosion q' vanishes in exactly 
the same way as the slope /i', or the profile h{x) is not 
differentiable at the crest. As the reader may already an- 
ticipate and will be verified below, both cases have their 
physical realizations, the former in heaps or small dunes 
with smooth crests and the latter in large dunes with a 
slip face that terminates in a sharp brink. 

The problem we face, if we want to calculate the dy- 
namic evolution of desert topographies, is the closure of 
Eq.(|l]) or (H) by expressing the flux q{x, t) in terms of the 
height profile h{x,t) and the external wind and bound- 
ary conditions. Since for the applications we have in 
mind, the migration velocity is very small compared to 
the speed of elementary sand transport processes (grain 
hopping etc.) and the wind speed, the topography can 
be assumed to be stationary for considerations concern- 
ing the wind and sand transport dynamics. This allows 
one to subdivide the problem of calculating q(x) into two 
independent steps. First, one needs to know the station- 
ary wind velocity above a given topography. More pre- 
cisely, what is required is the shear stress t exerted by 
the wind onto the ground. And secondly, one needs a 
model that predicts the stationary sand flux q{x) for a 
given stationary t{x), schematically 



h{x) 
t{x) 



t{x) 
q{x) 



(3) 
(4) 



Computing the derivative g' and integrating the mass 
conservation Eq. (^ then closes the model and allows one 
to predict the development of the surface profile in time. 
Since aeolian dunes typically have relatively gentle slopes 
outside their slip face, we will at the present stage restrict 
the scope of the minimal model to this case and disregard 
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in Eq.(^ the direct slope effects h'{x) — > q{x) onto the 
flux outside the slip face. 

In special cases, the relations (||), (jj) are phenomeno- 
logically and theoretically well established. For a flat 
surface, h{x) = const., it is well known ||2^ that the 
mean turbulent wind velocity increases logarithmically 
with height above the surface. It can be characterized 
by a single characteristic velocity, the "shear velocity" 
M, defined by ul = tq/ Qa with tq the (suitably time av- 
eraged) shear stress and Qa the density of air. Since the 
shear stress of the air is transmitted to the surface, the 
latter can mobilize grains on a surface covered with sand, 
if it exceeds a threshold value r*. As a result, the wind 
entrains some grains into a surface layer flow. The grains 
advance mainly by an irregular hopping process ("salta- 
tion"), thereby reducing the wind velocity in the surface 
layer. Via this feed-back mechanism a unique relation 
between the shear stress r and the sand flux q is estab- 
lished in the equilibrium state. If r is not too close to the 
threshold, this relation can approximately be represented 
as fl 



Qs oc r 



3/2 



(5) 



Although a host of more accurate descriptions have been 
discussed in the literature |l|, |l|, and one of 

them will be part of our definition of the minimal model 
below, the simpler Eq.(||) will be sufficient for our quali- 
tative discussion in the first part of the paper. The index 
s in Eq.(^ emphasizes that such local relations are re- 
stricted to situations where the flux is saturated, that is, 
equal to its equilibrium transport capacity. This is cer- 
tainly not the case near a boundary between uncovered 
and covered ground or on sloped beds. Neglecting this 
restriction for the moment, Eq.(^ predicts that the shear 
stress perturbation 



f (a;) = t{x)/tq ~ 1 



(6) 



above a modulated topography h{x) is responsible for 
flux gradients dqs/dx that cause erosion and deposition 
and thus — according to Eqs.(|l|), (||) — migration of the 
sand surface. Explicitly closing the model by assuming 
that the shear stress is an affine function of the modula- 
tion of the topography (f cx h) leads to what we call the 
"zeroth order" model, which will briefly be analyzed in 
the next paragraph. 



B. The "zeroth order" model 



The zeroth order model is given by 



f{h{x)} 
q{T{x)} 



f{h) (X h{x)/L 
q{T) = qs{T) 



(7) 
(8) 



where we have used the curly brackets to indicate a gen- 
eral functional dependence and introduced a character- 
istic length scale L of the topography to normalize the 



height profile. (The motivation for the latter step will be- 
come clear in the next section.) The zeroth order model 
assumes local relations in Eqs.(|^),(^. It approximates 
the wind shear stress perturbation by its "affine" con- 
tribution (proportional to the profile h that causes the 
perturbation) and replaces the true sand flux q by its 
saturated value g^, thereby neglecting saturation tran- 
sients. This model is so simple that its qualitative pre- 
dictions for an arbitrary smooth heap of sand can easily 
be anticipated without doing any actual calculations. 

Combining Eq.(|) with Eqs.(|)-(|) one obtains a sur- 
face velocity that increases with height {dv/dh > 0) 
due to the nonlinearity of Eq.(||). This implies that 
the upwind (or "stoss") slope tends to decrease and 
the downwind (or "lee") slope tends to increase. Since 
dq/dx cx dh/dx by the chain rule, there is no erosion or 
deposition at the top of a smooth heap, which therefore 
keeps its initial height. Obviously, integrating forward 
in time will eventually increase the lee slope up to the 
angle of repose, where surface avalanches have to be in- 
troduced and a slip face of constant slope develops. If the 
latter reaches the crest the above argument for the per- 
sistence of the height can no longer be applied, because 
the slope at the crest is then ill defined. Since there is 
so far nothing to stop a further decrease of the wind- 
ward slope, the model dune will then start to decrease in 
height and finally flatten out. The steady-state solution 
is a flat surface. 



The simple argument shows that the zeroth order 
model — although it gives some clue as to the origin of 
the slip face — is insufficient for a proper qualitative un- 
derstanding of dunes. However, some important lessons 
can be learned from it that will be helpful in our further 
investigation of the problem. First, even with a very 
simplistic model any reasonably heap-like initial condi- 
tion will quickly develop into a dune-like shape with a 
slip face. Secondly, although the latter may seem to con- 
verge to a steady-state solution for intermediate times, it 
finally turns out to be unstable and fiattens out. The dis- 
cussion of the migration velocity in Section II A suggests 
that small deviations from Eqs.(0) and (^) at the brink 
can make an important difference. Obviously some cau- 
tion is needed in judging the success of numerical models 
of dune formation. Unless stability has explicitly been 
demonstrated, they may be suspected to fail in a similar 
way as the zeroth order model when integrated over suf- 
ficiently long times (which has actually not been checked 
for some models that can be found in the literature) or 
to be sensitive to numerical errors at the brink. Detailed 
numerical modeling should therefore be preceded by a 
sound qualitative understanding of the mechanisms un- 



derlying dune formation. We will argue in Sections III 



[V that to this end a subtle balance between two small 
deviations from Eqs.(0),(^) and especially non-local con- 
tributions in Eq. (0) have to be taken into account. 
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C. Migration velocity 

Before entering a detailed discussion of the minimal 
model, it is worth pausing for some general thoughts as 
to what can be said about the shear stress and the speed- 
up of the wind above an obstacle, without actually doing 
the (somewhat involved) calculation. 

A basic property of strongly developed turbulence is 
its dilation invariance or scale free structure. Whereas 
general Navier-Stokes flow is invariant under a scale 
transformation that keeps the Reynolds number con- 
stant, strongly turbulent flow (for "infinite" Reynolds 
number) allows for infinitely many such similarity trans- 
formations. Landau and Lifshitz took advantage of 
this fact for deriving the logarithmic velocity profile men- 
tioned above by an elegant scaling argument. The log- 
arithmic velocity profile suggests that the speed-up of 
the wind and therefore also the shear stress perturbation 
above a heap of given shape should itself be logarithmi- 
cally dependent on its size. But how do they depend 
on the shape of the obstacle? Since the flow itself does 
not provide any characteristic length scale, the dimen- 
sionless quantity f defined in Eq.(^) can only depend on 
a dimensionless characterization of the profile h{x). In 
other words, to lowest order in the perturbation, it must 
be a linear functional of the derivative h' and can be 
written as 

f(C)=er{/'(0} , e^H/L, (9) 

with a dimensionless profile function 

fiO ^ h{x)/H C^x/L. (10) 

and a scale-free (and necessarily non-local) linear func- 
tional T. This reasoning can be repeated for the dimen- 
sionless velocity and pressure perturbations. Intuitively, 
the scaling f cx e for a flat smooth obstacle (e ^ 1) can 
be understood from Fig. |^. When the air flows over the 
obstacle, the velocity close to the obstacle is deflected by 
an angle e whereas it remains constant far above the ob- 
stacle. For incompressible flow, continuity translates this 
into a speed-up of order e and (via Bernoulli's law) into 
a corresponding pressure drop near the top of the heap. 
This in turn causes a shear stress perturbation f of the 
same order. 

These general considerations already allow us to pre- 
dict the scaling of the migration velocity v with dune size 
if we assume that dunes of different size have roughly 
similar shapes /(^) and aspect ratios e, which is indeed 
suggested by the scale invariance of the turbulent wind 
field and by observations. Inserting Eq.(^) into Eq.(||), 
and again approximating Eq. (^) by a local sand transport 
law q = qij), we find 

dq dr/dx ^ dq T' {f} 1_ 
"df g.ef' df QsLf L- ^ ' 

The final proportionality strictly holds only if the steady- 
state shape /(^) and aspect ratio e are scale invari- 
ant. However, it can be expected to be robust and 
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FIG. 2: Schematic sketch of the deflection of the wind veloc- 
ity u above a flat heap of aspect ratio e = H/L <^ 1. The 
characteristic length scale L is in this context conventionally 
often identified with the half length at half height of the heap. 
The vertical deflection causes a speed-up above the top of the 
heap. This is accompanied by a pressure perturbation that is 
negative above the top of the heap and positive at its tails. 
Due to turbulence, the flow pattern is asymmetric even above 
a symmetric heap. 

rather insensitive against violations of exact scale invari- 
ance. First, the normalized steady-state shapes /(^) are 
strongly constrained by the requirement that they render 
v{x) = V, independent of x, along the heap. Therefore, 
they should to a first approximation be independent of 
size, which is indeed borne out by the minimal model 
(Fig. |l^) and empirical observations Moreover, the 
dependence v{f'{^)} is rather indirect and can therefore 
be expected to be weak. Secondly, Eq.(^) suggests that 
for gently sloped obstacles (e <C 1) the dependence of 
dq/dr on the aspect ratio s also is not very pronounced. 
And finally, — due to the above mentioned scale invari- 
ance of turbulence — a scale invariant aspect ratio can 
reasonably be expected for large dunes. In fact, we will 
show below that the minimal model predicts that the 
aspect ratio of small heaps is not constant but rather de- 
creases proportional to their height. But this also implies 
that the latter becomes too small to have a very signifi- 
cant effect on the above argument. Note, however, that 
only for strictly scale invariant dunes, Eq. (pT|) becomes 
identical to the often quoted observation that dunes mi- 
grate with a speed inversely proportional to their height 
pl[ . Since the deviations of large dunes from scale in- 
variance are not very pronounced, the difference between 
these predictions is not very strong except for small dunes 
and heaps. Presently available field data are maybe not 
accurate enough to clearly distinguish between v oc 1/L 
and V (X 1/H, though some data support v cx 1/L, most 
notably the comprehensive study of barchan dunes in 
southern Peru by Finkel |2^. As we will show below, 
our numerical results for the minimal model clearly favor 
V oc 1/L. 

We also note that together with Eq.(^), Eq.(|l^) more- 
over predicts that the migration velocity grows non- 
linearly with (as the third power of) the wind velocity. A 
more accurate relation can be obtained from the minimal 
model as described below, but the qualitative conclusion 
is the same. Dunes can migrate farther in a short period 
of exceptionally strong wind than during much longer 
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periods of gentle winds. Finally, we should mention that 
some caution is needed when identifying the characteris- 
tic length scale L in Eq.(|l]). In our discussion, we have 
so far assumed that /(^) is a smooth function, which 
is not the case for dunes with slip face. Below we will 
argue that in this case / should be identified with the 
envelope of the dune and its separation bubble and L 
with the characteristic length scale of this envelope. For 
a barchan dune the latter practically coincides with the 
total length of the dune from its windward end to the 
tips of its horns (cf. Fig. . 

In contrast to the overall migration velocity of a trans- 
lation invariant dune, the position dependent migration 
velocity v(x) that determines that shape is much harder 
to obtain since it requires a precise knowledge of the non- 
local functional T in Eq.(^). This will be provided in the 
following section. 
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FIG. 3: The theoretical prediction for the dependence of the 
parameters A (solid line) and B (dashed line) of Eq.(p^) on 
the ratio of the dune size L to the roughness length zq. 



III. WIND SHEAR STRESS 

A. Surface shear stress on a smooth heap 

The discussion in the preceding paragraph showed that 
— in contrast to the assumption in Eq. (|^) — the depen- 
dence of the shear stress on the height profile is non- 
local. Although it will turn out that this shortcoming of 
Eq.(0) is not responsible for the failure of the zeroth or- 
der model, it should by now have become apparent that 
further progress can hardly be achieved without a rather 
detailed understanding of the turbulent wind field above 
heaps and dunes. For dunes with a slip face that typically 
has a slope of about 32° — 35° and terminates in a sharp 
brink, the situation is similar to the textbook example 
of a backward facing step, which has the reputation of a 
test case for numerical turbulent models. Even if a com- 
mercial turbulent solver is used, the accurate calculation 
of the shear stress e.g. on a barchan dune is a non-trivial 
task and quite demanding in computer time and memory, 
and the most interesting long-time dynamics of dunes is 
therefore difficult to access. For this reason, we want 
to focus on flat smooth heaps, first. In this case, one 
can apply an analytical perturbation theory for turbu- 
lent boundary layer flow over smooth hills that has been 
developed over the last decades |l[ [l2[ 0. 
Though the calculation is essentially a formalization of 
the intuitive description accompanying Fig. ^ it requires 
a highly non-trivial boundary layer construction that we 
will not recapitulate here. The interested reader is re- 
ferred to the original literature. We merely quote the fi- 
nal result for the a;— component (along the main wind di- 
rection) Tx of the surface shear stress perturbation above 
a profile h{x,y) 

^^y{r^} = ^M-, y)} ■ (12) 



(fc2+fc2)l/2 



ky by Txy- For simplicity the logarithmic A:— dependence 
of the parameters A and B was neglected. The latter are 
then given by 



ln($2/ln$) 
^ " 2(ln0)3 [1 + + 2 ln(^/2) + 47] 

B tt/ [1 + In + 2 ln(7r/2) + 47^] 
(j)= 2k^^/\tl4) 



(13) 



We have abbreviated the Fourier transformation from the 
space variables x, y to the respective wave numbers k^^ 



and depend logarithmically on the ratio $ = L/ zq, where 
L is the characteristic length of h{x, 0) (for this pur- 
pose conventionally often identified with half the length 
at half height or about one fourth of the characteristic 
wavelength) and Zq is a measure of the surface rough- 
ness (typically an effective length somewhat below the 
linear dimension of the latter). We also have introduced 
the von Karman constant k, k, 0.4 and Euler's constant 
~ 0.577. A practical approximation for (j) is obtained 
by iterating (twice) the implicit equation for (j) and clos- 
ing it by dropping the remaining ln0. The dependence 
of A and i? on L is depicted in Fig. ||. Obviously, as 
long as L/zq does not change by orders of magnitude 
(e.g. due to vegetation), this extremely weak scale de- 
pendence is negligible for our purposes, and A and B 
can be regarded as constant theoretical or phenomeno- 
logical parameters. For dcfinitcness we often work with 
the values A — A and B — 0.25 (approximately obtained 
for L/zq = 10^ . . . 10^). Although these values may dif- 
fer somewhat depending on the particular application in 
mind, or on the presence or absence of ripples, and may 
phenomenologically be somewhat different from the the- 
oretical prediction, this docs not affect our general con- 
clusions. 

For the following discussion we want to specialize 
Eq.(|l2[) to the central slice of a transverse or barchan 
dune along the wind direction. To this end we evaluate 
Eq.(p^ for the central slice h{x) of a heap that has a 
Gaussian shape with standard deviation a in the trans- 
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verse direction parameterized by y, 

h{x,y) = h{x)e-y'/'^''^ 



(14) 



This approximation is technically useful, and although it 
may seem relatively crude for a particular real dune, it 
typically does not introduce any noteworthy derogation 
of the results compared to a more accurate description. 
The Fourier coefficients of the shear stress f{x) = fx{x, 0) 
on the central slice along the wind direction become 



a A 

T{t{x)} = -^fc(fc + iB\k\) X 



e 4 



T{hix)} (15) 



Here, Kq denotes a modified Bessel function and the 
Fourier transforms are one-dimensional, so that we can 
drop the redundant x— subscripts. 

For transverse dunes (cr/L — s- oo), we obtain the two 
equivalent expressions 

T{t°°{x)} = A{\k\ + iBk)T{h{x)} . (16a) 
f°°(x) = A [h'{x) (g) (ttx)-^ + Bh'{x)] . (16b) 

For the real space version we have abbreviated a convo- 
lution integral according to 



(17) 



Evaluation for arbitrary a gives two correction terms 

f'^ =f°°- A(/i® Ai+B/i'® A2) , (18) 

with 



U 



Ai 
A2 



(3 ^ ^u(^ 



1 



— / a4 cos — 
an Jo \a 



TTX^ 



(19) 
(20) 



two even functions depicted in Fig. ^ that are flat for 
a/ L — !■ 00 and become peaked for a L. (The confluent 
hypergeometric U functions have been introduced to 
rephrase the sine-part of the Fourier integrals.) 

Since the correction terms in Eg. ([l^) are numerically 
small, we may — given a reasonable localized heap shape 
in the wind direction — approximately replace both func- 
tions Ai and A2 by delta functions, thus arriving at 



Too - A[ci{(j)h + Bc2{a)h'] 



(21) 



In this approximation they are seen to give merely a 
(T— dependent renormalization of the asymmetry param- 
eter B B{a) < B{oo) — B and to add a ( trivial) 
term ci{a)h{x) within the brackets of Eq.( |l6b| ). Nu- 
merically, one can estimate Lci{L/^/2) and C2{L/^/2) 



Ai 




. 05 



A2 




FIG. 4: The peaked functions Ai and A2 of Eq.(na) for cr = 1, 
2, 5. The area under the peaks remains constant (0 and l/27r), 
while the peak heights decrease proportional to cr"^ and , 
respectively. 



to be about 0.2 (cf. Fig. ||). The exact cr— dependence 
of the coefficients is determined by the shape and ex- 
tension of the heap in the a;— direction, because the area 
under the peaks Ai and A2 is constant and independent 
of (7, while the peak height decreases proportional to <j^^ 
and (7^^, respectively. Since both corrections vanish for 
a/ L 00 and do not contribute any substantial new 
effects to Eq.(p^, they may for simplicity be omitted al- 
together in the following discussion that mainly aims at 
a qualitative understanding. Fig. |^ moreover shows that 
they can approximately be mimicked by a renormalized 
parameter A in Eq. (|l^) for the central slice of a symmet- 
ric heap. This leads to the important conclusion that the 
wind shear stress on the central slice of a 3d symmetric 
heap and on a heap with a profile that is constant in 
the transverse direction, are qualitatively the same and 
quantitatively similar, which was not a priori obvious. 
Together with the fact that on a gently sloped obstacle 
the transverse components of the shear stress are small 
compared to its longitudinal components, this suggests 
that the predictions of Eq. (|l6|) apply in a first approxi- 
mation to any slice of a dune parallel to the wind direc- 
tion. In this sense, the study of Eq.(p^ is representative. 
Summarizing the foregoing discussion we can say that to 
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FIG. 5: Shear stress perturbation above the central slice of a 
Sd symmetric (a = L/\/2) Gaussian heap. The plot compares 
two approximations to Eqs.(^, ( |l8| ) (points): (i) Eq.(|2l|) 
with Lci « C2 ~ 0.2 (solid line), and (ii) Eq.([La) with A renor- 
malized by a factor 0.8 (dashed line). While (i) is practically 
indistinguishable of Eq.(^) on the present level of accuracy, 
the simpler approximation (ii) already captures most of the 
3d effects. 



gain a qualitative understanding of dune formation by 
aeolian sand transport one may focus on Eq.(|l^) as a 
model for the wind shear stress. We therefore analyze 
this equation in some detail in the next paragraph. 



B. Properties and consequences of Eq.(16) 



A scaling analysis of Eq. (^^ immediately reveals that 
f is indeed of the general form anticipated on general 
grounds in Eq.(H). The amplification of the shear stress 
at the top of a smooth profile is thus determined by its 
aspect ratio e = H/L and is essentially independent of 
the absolute height H . It only has a very weak log- 
arithmic dependence on the absolute size of the dune 
through the prefactors A and B given in Eq. (|l^) . More- 
over, for a symmetric profile /{—£,) = /(C): ^he sum 
of a symmetric part and an anti-symmetric part, i.e., 
the flow over the heap has a symmetry breaking com- 
ponent that is a consequence of turbulence. The origin 
of the symmetric and antisymmetric parts of f can intu- 
itively be understood as follows. As we have pointed out 
in Section II C (see Fig. ||), the streamlines have to be 
compressed above the heap if the perturbation is not to 
be transmitted to infinite height, and as a consequence, 
there is a corresponding increase in the shear stress. For 
the laminar average flow, this speed-up and the associ- 
ated decrease in atmospheric pressure above the heap are 
symmetric for a symmetric heap as is the corresponding 
shear stress perturbation, which accounts for the domi- 
nant symmetric part of f. On the other hand, the inertia 
of the turbulent velocity fluctuations around this lami- 
nar main flow contributes an asymmetric resistance to 
deflections of the flow. It counteracts the upturn of the 



FIG. 6: Lower curves: The normalized profiles of Eq.(p2[), 
the Lorentzian /l (dashed), Gaussian fc (solid), and cosine 
(dotted). Upper curves: The corresponding surface shear 
stresses t{x)/to from Eq.Q with Ae = 0.8, B = 0.25. 



streamlines on the windward side and the downturn on 
the lee side. Formally, this effect enters the perturbative 
calculation of f through the Reynolds stress. 

Further insight can be gained from special analytical 
solutions to Eq.(p^. For the normalized heap profiles 



/G(C)=exp(-C2) 
/S(0 = ^(C)cos"C 



(22) 



with 



[l m<7r/2 
10 m>n/2 



(23) 



we obtain 



TG = 2A 
A 



-1/2 



-e(S + erfiO/G(0 



fc = - cos(20 [ Si(^ + 20 + Si(7r - 2^] ^^4) 

- A sin(20 [ Ci(^ + 20 + Ci(-^ - 20 

-Ci(7r - 2a;) - Ci(-7r -f 22;)] 
-2BS'(0 cosC sin^ (n = 2) 

with erfi the imaginary error function and Si and Ci the 
sine and cosine integral functions, respectively p3t . The 
result given for fc is for the special case n = 2. Both 
the pro files of Eq. (|2^) and the corresponding solutions of 
Eq.(|l6[) given in Eq. (|24|) are shown in Fig. |g. 

The plots of T show that as a rule of thumb one can 
estimate the relative magnitude of the shear stress per- 
turbation at the top of the heap by Ae. The plots share 
several crucial properties not present in the affine approx- 
imation f cx ft, of the zeroth order model. At the tails of 
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FIG. 7; The symmetric and asymmetric parts fsym and fasy of 
the shear stress perturbation tg of Eq.(2^ for the Gaussian 
profile /g of Eq.Q with As = 0.8 and B = 0.25. Note 
the small windwara shift of the maximum of the shear stress 
with respect to the crest of the heap caused by the asymmetric 
contribution proportional to B. 



the profiles in Fig. |[ the shear stress decreases below its 
asymptotic value tq on the plane. This effect is partic- 
ularly pronounced for the profile that has a disconti- 
nuity in its second derivative. Further, as a consequence 
of the second term in Eq. (^6|) , the surface shear stress is 
not symmetric even for symmetric profiles like those in 
Eq. (p2[) . Fig. displays the symmetric and asymmetric 
parts fgym and fasy of the shear stress perturbation tg 
for the Gaussian profile fa, separately. The asymmet- 
ric contribution to f is small compared to the symmetric 
one. For the profile Jl the corresponding shift Sxr of 
the location of the maximum of the shear stress with re- 
spect to that of the maximum of fL{x) can be calculated 
analytically, 



Sxr/L = 2 (1 + ^2)1/2 sin[arctan(B)/3] - B . (25) 

It is indeed found to be very small, because B is small 
and thus Sx-r/L « —B/i typically amounts to a length 
of about a few percent of the total heap length. Never- 
theless, it is a crucial element in the modeling of aeolian 
sand transport, as will now be demonstrated. 

For a qualitative estimate of the effects of Eq.(p^) 
onto the sand transport over a dune, it is useful to con- 
sider once again the local zeroth order model for aeo- 
lian sand transport, Eq.(^), i.e. a completely saturated 
flux q = (/s(t) with given by Eq.(|^). (Below, we will 
show that this is asymptotically valid on large dunes in 
strong winds.) The distinct features of Eq.([l6[) that are 
missing in the zeroth order approximation for the shear 
stress, Eq.(|^), are then easily seen to have potentially 
profound effects on the shape evolution. First, due to the 
depression of the shear stress at the tails of the profiles, 
deposition rather than erosion may occur at the wind- 
ward foot. Secondly, due to the asymmetric contribution 
in t{x) there can be a net deposition on a symmetric 



v{x) 




FIG. 8: The position dependent surface migration velocity 
v{x) in arbitrary units according to Eqs.(|2|), ( |l6| ) with A = A, 
B = 0.25. Upper part: For the profile /q (gray) and varying 
aspect ratios e = 0.01 (dashed), 0.1 (solid), 0.19 (dotted) and 
the local flux relation Eq.(H). Lower part: For the profile /^''^ 
(gray) with e = 0.1 (solid), 0.05 (dashed) and the non-local 
flux relation Eq.(po|). For simplicity, qs was represented by 
Eq.(|^ and ~ 0.1 was taken constant. 



heap of sand. In particular, the shift of the position of 
the maximum shear stress with respect to the top of the 
heap allows deposition at the top of the heap. For an 
initially flat heap of sand there is thus the possibility 
of a steepening of the windward slope and mass growth. 
This implies that a plane sand surface is unstable against 
modulations. To illustrate this effect, we used Eq.(p4|) to 
calculate the migration velocity v{x) of a cosine-shaped 
heap of sand f^ix) according to Eqs.(||), (||), (||), and 
(|l6|). The latter is shown in the upper part of Fig. ||. 
The decrease of v{x) on the lee side reveals the antici- 
pated self-amplifying tendency of the unstable lee slope 
to steepen, since v' < implies that v increases with 
height. This gives rise to the formation of the slip face. 
More interestingly, Eq.(p^ renders v{x) approximately 
constant over almost the whole windward side if £ = H/L 
is close to a certain value determined by the values of 
the coefficients A and B in Eq.(p^. Slightly better con- 
stancy can be achieved for slightly lower n (with slightly 
larger e) but not for the profiles fa and /l, for which 
v{x) is always non-uniform. Let us finally consider the 
dashed and dotted lines in the upper part of Fig. ||. They 
were obtained for a smaller and a larger aspect ratio and 
represent a steepening {v' < 0) and flattening {v' > 0) of 
the windward side, respectively. In other words, profiles 
with larger/smaller windward slopes are driven towards 
the solution with constant windward v{x) and a stable 
optimum windward slope different from zero. Altogether, 
Fig. H thus suggests that the coupled Eqs.(||) and ( [l^ ) 
drive a heap of sand towards a "dune" with a cosine-like 



windward profile of a preferred aspect ratio, and a slip 
face on the lee side. 

From the qualitative analysis presented so far, it is not 
yet obvious that this process converges to a translation 
invariant steady-state solution. Several previous studies 
using similar descriptions either did not scrutinize the 
long time behavior of their models or failed to 

obtain stable dunes |]l3[ ^ . To obtain a consistent gen- 
eral model for dune formation under general influx and 
wind conditions, the present wind model Eqs.(p^, ( p^ ) 
has to be appropriately adapted to situations with flow 
separation above slip faces. And, most importantly, the 
saturated-flux approximation Eq.(H) of the "zeroth or- 
der" model has to be abandoned. These steps will be 



streamlines 



separating streamline 

separation bubble 



discussed in the following subsection and in Section IV, 
where we will also explain the lower part of Fig ^. This 
will complete the definition of the minimal model. Its 
numerical solutions will be presented in Section 



C. Flow separation 

The wind model as discussed so far works fine for 
smooth heaps with gentle slopes. However, as we have al- 
ready mentioned, its application to dune profiles with slip 
faces and sharp brink lines is not straightforward. The 
perturbative turbulent boundary layer approach leading 
to Eq.(p^ does not account for flow separation, a phe- 
nomenon that occurs at sharp edges and steep slopes 
to prevent an extreme bending of the streamlines [ pO[ . 
(For some of the technical terms involved in this sec- 
tion, the reader is referred to Fig. |[) Instead of bending 
the streamlines around sharp edges, re-circulating ed- 
dies separate from the (on average) laminar main flow, 
thereby creating an effective envelope that diverts the 
main flow on a smooth detour around the obstacle. See 
Figs. H and |l^ for a schematic sketch and a numerical 
calculation of a typical velocity fleld, respectively. For- 
tunately, it turns out that dune formation and migration 
do not in general depend very sensitively on the details 
of this complicated process. Or in other words, there is 
a large number of interesting problems of aeolian sand 
transport for which these details are largely irrelevant, 
and for which their somewhat realistic physical represen- 
tation would create a huge overhead in complexity (es- 
pecially in 3c?) to an otherwise tractable problem. It was 
therefore suggested earlier that for the purpose of cal- 
culating the shear stress on the windward side of a dune, 
one may to a good approximation represent flow separa- 
tion on the lee side by the following heuristic method. A 
wind model such as Eq.(p^) restricted to smooth, gently 
sloped objects is applied to the envelope 



h{x) 



i{h{x),s{x)} 



(26) 




of a dune h(x) and a phenomenologically deflned separa- 
tion bubble s(x). This disregards the fact that the sepa- 
rating streamline does not represent a solid boundary of 



FIG. 9: Sketch of the central slice of a barchan dune and the 
separation bubble. The shear stress on the windward side of 
the dune is calculated by applying Eq.(jl^) to the phenomeno- 
logically defined envelope of the dune and the separation zone. 



the same roughness as the original object, but the corre- 
sponding errors are expected to be small. Typically, one 
wants s{x) to be a mathematically simple smooth con- 
tinuation of the dune proflle. It is, however, crucial that 
the latter respects some major phenomenological prop- 
erties of flow separation |2j. Although this is by no 
means a rigorous procedure, one can test its predictions 
for selected cases against numerical solutions of various 
turbulence models. The hope is that via this approach, 
one can eventually get a qualitative understanding of the 
mechanisms and phenomena involved in dune formation 
and migration, leaving certain quantitative aspects to a 
more elaborate (and much more laborious) future analy- 
sis. 

In the spirit of the minimal model we want to parame- 
terize the separating streamline s{x) in the simplest form 
that obeys physically motivated boundary conditions at 
its detachment and reattachment points Xd and Xr- At 
detachment, the slope of the separating streamline must 
match the slope of the dune. Moreover, also the curva- 
ture must be continuous there, since discontinuities in 
curvature are detected by Eq.(nq) and cause kinks in r 
and discontinuous steps in the erosion/deposition as is 
e.g. the case for the profile /q. For the reattachment 
point, there are no comparable restrictions to the slope 
and curvature, since the separation bubble is not sharply 
defined there, and the model aims at a realistic descrip- 
tion of the conditions in the wake region only insofar as 
they affect the shear stress on the windward side. On the 
lee side, inside the separation bubble, the shear stress 
can simply be set to zero [^8|, since it is typically be- 
low the threshold for aeolian sand transport. Therefore, 
the choice of the reattachment matching condition is a 
matter of convenience rather than physical significance 
in the present model. However, we want s{x) to repro- 
duce some common phenomenological knowledge about 
flow separation. First, from many numerical calculations 
it is known that, at high Reynolds numbers, the turbu- 
lent boundary layer reattaches at a distance of about 6H 
after a backward-facing step of height H. Secondly, it 
has often been observed experimentally that in strongly 
turbulent flows over hills and symmetric triangular ob- 
stacles, flow separation sets in if the backward slope ex- 
ceeds an angle of about 14°. Although, in both cases the 
exact numerical values depend on various factors such 
as the surface roughness and the Reynolds number, they 
shall be treated as fixed phenomenological constants at 
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FIG. 10: Separation bubbles with a maximum negative slope 
of 0.25 according to Eq.(p9[) for varying initial slopes —0.25 < 
h'd < 0.25. (The aspect ratio of the plot was stretched for 
presentation.) 



the present stage. A model that fulfills all the above re- 
quirements is a third order polynomial with continuous 
slopes at the boundaries and a maximum negative slope 
of tan 14°. The boundary conditions 



s{xd) 
s'ixd) 



h{xd) 
h'{xd) 



max{— s'(a::)} 



s{Xr) — 
s'{Xr) = . 

tan 14° = 0.25 



(27) 



constrain the third order bubble parameterization to be 
of the form 

s{z) = {2hd + h'^Lb)z^ - {3hd + 2h'^Lb)z^ + h'^L^z + hd 

(28) 

with z = {x — Xd)/ Li, e [0, 1]. With the further abbrevi- 
ation V = h'^/s\ we can express the length Lh = Xr — Xd 
of the bubble as 



3/irfl 



3hd 
2s' 



1 



(29) 



where the final approximation for small h[i is sufficient for 
our purpose (and numerically better behaved as the exact 
expression) . A subtlety of such a separation bubble pa- 
rameterization is the fact that the slope at Xd determines 
the length of the bubble, which in turn, via Eq.(p^) influ- 
ences the curvature at Xd- In other words, the presence of 
the bubble introduces a non-local feed-back between the 
slope and the curvature at the brink, which we believe is 
physically reasonable. In Fig. |l^ we give some examples 
of separation bubbles for different boundary slopes /i^, 
while Fig. |ll| illustrates the application of the above dis- 
cussion for the calculation of the shear stress. It shows an 
example for a dune profile h{x) with slip face and the sep- 
aration bubble s(x), together with the shear stress t{x) 
resulting from Eq.(p^ if h is replaced by the envelope h. 

We have performed several series of numerical fluid 
dynamics calculations in 2c? and 3d with the commercial 




x/L 

FIG. 11: The windward profile h{x) of a dune with slip face 
and the separation bubble s{x) form together a smooth ef- 
fective obstacle, defined by the envelope h{x). To calculate 
the shear stress t{x) on the windward side of the dune, h is 
substituted for h in Eq.(|l^). In the region of re-circulation 
the surface shear stress r is set to zero |28| . Without the sep- 
aration bubble, t{x) would develop a sharp singularity at the 
brink. 



fluid dynamics solver Fluent 5 ||29| using the fee and large- 
eddy turbulent closure models to conflrm the general pic- 
ture outlined above and our particular implementation of 
the separation bubble. The differences between numeri- 
cal and theoretical predictions for the shear stress on the 
windward side of various dune- and heap-like objects in 
2d and 3d were quantitatively small and not more signif- 
icant than other neglected terms. Moreover, a compari- 
son of predictions obtained from Eq. (|2^) with wind mea- 
surements on a barchan dune in Brazil pO| showed good 
agreement. Therefore, we are confldent that the pro- 
posed mathematical description of the wind shear stress 
captures the relevant aspects in the spirit of the mini- 
mal model. As an example for the numerical fluid dy- 
namics calculation, we show in Fig. ^ the flow velocity 
in the symmetry plane of a 3c? barchan dune obtained 
with Fluent 5 The wind is blowing from left to 

right. The boundaries were chosen to be periodic in the 
transverse direction. At the influx boundary, the veloc- 
ity was fixed by imposing a logarithmic velocity profile. 
The wind profile at the outflux boundary is not known 
a priori. Although, for high Reynolds numbers the lat- 
ter is expected to affect the solution only close to the 
boundary, it is well known that different choices for the 
outflux boundary condition as well as different discretiza- 
tion schemes may lead to quantitatively different results 
Here, we chose to set the derivative of the veloc- 
ity normal to the outflux boundary to zero. The surface 
profile was represented as a solid boundary with constant 
roughness length. Finally, along the top boundary we im- 
posed the velocity of the undisturbed logarithmic inflow 
profile. The whole calculation was performed on a grid 
that had an exponentially growing mesh size in the ver- 
tical direction. A considerable grid refinement was nec- 
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FIG. 12: Cut along the symmetry plane of a 3d barchan dune. 
The velocity vectors calculated numerically with a commercial 
fluid dynamics solver clearly display the flow separation 
at the brink and a large eddy in the wake region. 



essary in the wedge-like region of the separation bubble 
close to the brink. 

These remarks complete the first task of constructing 
a model for the calculation of the wind shear stress on 
a given dune profile as outlined in Eq.(H). By deriving 
the linear Eq.(^6|) for the shear stress and combining it 
with the heuristic separation bubble, we have obtained an 
approximate but numerically extremely efficient model 
for the wind shear stress on dunes. This is a crucial step 
in the construction of a minimal model of aeolian sand 
dunes, since the enormous complexity of the turbulent air 
flow over structured terrain otherwise severely restricts 
the possible applications of the model. 

Going back to the upper part of Fig. || with the above 
discussion in mind, we can re-interpret this figure in or- 
der to anticipate the behavior of the surface migration 
velocity v{x) of a dune with slip face. If, for qualita- 
tive purposes, is interpreted as the envelope of a dune 
and its separation bubble, we can conclude that the slip 
face must be located near the sharp drop of v{x) slightly 
upwind from the top of the envelope. This is indeed 
consistent with observations for large dunes. Together 
with the good representation of the windward profiles of 
large dunes by (n w 2), it suggests that the given 
description becomes qualitatively correct in the limit of 
large dune sizes. The next section is devoted to the dis- 
cussion of important subtleties related to the fact that 
dunes are not typically in this limit. 



IV. SAND FLUX 

As outlined in Eq.(^, the second task in the speci- 
fication of the minimal model is to find a prescription 
for calculating the sand flux q(x) for a given topogra- 
phy h{x) and shear stress t{x). So far, we have been 
using the local saturated-flux approximation Eq.(||) in 
our qualitative arguments. However, a closer look at the 
predictions obtained within this approximation reveals 
a number of inconsistencies. First, as we have already 
noted in the discussion of Fig. ||, the use of EqiH) to- 
gether with the complete wind model of Section |ll leads 
to the odd prediction of deposition at the windward foot 



of an isolated heap or dune, where the shear stress de- 
creases. This defect of Eq.(||) has been noticed in the 
literature before (see e.g. Refs. [^0[ Q). Previous nu- 
merical studies tried to avoid this problem by focusing 
onto the short time behavior and by introducing ad hoc 
heuristic methods such as a "smoothing operator" or 
an "adaptation length" . The reason for the problem 
is that the saturated-flux approximation breaks down at 
a boundary ground/sand. As another shortcoming, we 
want to mention that the model as discussed so far pre- 
dicts a universal scale invariant dune shape with a brink 
that is displaced slightly upwind from the maximum of 
the envelope, leading always to a positive slope at the 
brink. A glance at a real dune field proves that the 
latter is not always the case and careful measurements 
have revealed systematic deviations from scale in- 
variance. Though less obvious, it turns out that the rea- 
son for this discrepancy lies again in the saturated-flux 
approximation. Both mentioned problems are thus natu- 
rally resolved by introducing a slightly more general sand 
transport law that allows for saturation transients. 



A. Saturation transients 

The saturated-flux approximation Eq. (||) assumes that 
the flux is everywhere equal to the equilibrium transport 
capacity Qs of the wind. However, due to variable wind 
or sand conditions, the actual sand flux q is in general 
different from Qs- These deviations are called saturation 
transients, because they quickly relax to zero under ho- 
mogeneous conditions. We have recently demonstrated 
|]lO| that this relaxation occurs within a characteristic 
length scale, called the saturation length i^, which is re- 
lated to (but distinct from) the mean saltation length 
of the grains. It was moreover shown how the introduc- 
tion of saturation transients cures the problem of depo- 
sition at the windward foot of an isolated sand dune. 
Here, we only summarize the most pertinent results of 
this earlier development in order to demonstrate how a 
size dependence of the dune shape naturally results as a 
consequence of saturation transients. 

The sand transport model of Ref Q is based on a 
mean-field like description of saltation. It models a typ- 
ical grain that is accelerated by friction with the air and 
slowed down by dissipative interactions with the bed. 
The average properties of the complicated splash process 
p^ , ^ ^ are subsumed into two dimensionless param- 
eters, an effective restitution coefficient a for collisions 
with the bed, and a kinetic coefficient 7 that character- 
izes the relaxation of the density of saltating grains to 
its saturated value. Together with an effective height for 
the wind-grain interaction that enters only logarithmi- 
cally, these are the only phenomenological parameters of 
the model. They have been determined by a comparison 
with experiments and grain scale simulations. Formally, 
the model consists of two coupled differential equations 
for mass and momentum conservation, and a modified 
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turbulent closure relation that accounts for the feedback 
of the saltating grains on the wind velocity. 

For the present purpose, the model can be simplified 
by taking advantage of the fact that the prevailing condi- 
tions in applications to dunes are typically well described 
by the steady-state {d/dt ~ 0) version. Further, the re- 
laxation of the typical sand transport velocity can be 
assumed to be fast compared to the variations in the 
density of mobilized grains in the saltation layer. Ap- 
proximating the latter by its saturated value for the cal- 
culation of the effective wind speed Wog via the modified 
turbulent closure, one can decouple the mass and mo- 
mentum conservation equations. The whole model can 
then in a reasonable approximation be reduced to a sin- 
gle differential equation 

Isdq/dx = q{\ - q/q,) (30) 

for the sand flux q{x). The shear stress dependent pa- 
rameters 
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FIG. 13: The saturation length Is in meters as a function of 
the shear stress exerted by the wind onto the sand bed. This 
function sets the natural length scale for dunes and heaps. 



Is = l/ir/Tt - 1) 



PsUs 



(31) 



are immediately identified as the saturation length and 
the saturated flux, respectively. The equation for q^ gen- 
eralizes Eq.(|^) to arbitrary wind speeds. In the following 
we specify the explicit expressions for both quantities as 
they result from the sand transport model of Ref. [p^ , 
but the structure of Eq. ( pO[ ) is thought to be more general 
and independent of the precise form of Eq.(|3l|). Again, 
t(x) is the position dependent shear stress discussed in 
Section III and tj « 0.1 kgm^^s^^ is the estimated im- 
pact shear stress threshold that corresponds to a critical 
shear velocity u^t ~ 0.28 ms~^ (For simplicity, we 

do not introduce the additional threshold for purely aero- 
dynamic entrainment here, but allow instead for a small 
residual influx even if the latter is nominally zero.) To 
make the underlying structure of the model more palpa- 
ble, we have expressed Is in terms of another character- 
istic length scale / = 2au^/((77), which (up to a numer- 
ical factor) is the average saltation length of the grains. 
The latter — but not ts — must always be considerably 
smaller than the dune length for the model to be applica- 
ble. Further, we have decomposed qs into the saturated 
density ps = 2a(r — rt)/f; and the effective sand trans- 
port velocity at saturation Ug = Weff — <5w with UoS the 
effective wind velocity that accelerates the grains, given 
by 



- = 2yjTt + (r - rt)/C + (InC' - 2) 



(32) 



By g we have denoted the gravitational acceleration and 
from Ref ]lO| we adopt the (approximate) numerical val- 
ues a = 0.35, 7 = 0.2, ( = 8, C' = 200, and 8u = 1.8 m/s 
for the lag velocity of the grains. We note that these 
numerical values are not completely independent of each 
other and of the mentioned value for the impact thresh- 
old Tj, due to the calibration of the sand transport model 
with experimental data pO|. For convenience we show a 



plot of the saturation length Ig obtained with these val- 
ues in Fig. |l^. This completes the definition of the sand 
transport part of the minimal model on gently sloped 
ground. 



B. Consequences 

Before we complete the general model definition by 
a brief paragraph on slip faces, we want to point out 
some implications of the model as developed so far. First, 
note that the full expression for qs given in Eq. (|3l|) , con- 
tains Eq.(|^) as a limiting case for strong winds but is 
better approximated by cx t — tj for moderate wind 
speeds. For weak flux gradients and strong winds, one 
may set the left hand side of Eq.(pO|) to zero, leading to 
q — qs- This is typically the case on most of the windward 
slope of a large dune, where the left hand side of Eq.(|30|) 
can roughly be estimated by (sQs/L ^ qs- The local 
saturated-flux approximation with Eq.(^ for qs, which 
we have applied throughout our qualitative discussion 
so far, is thus asymptotically valid for large dunes and 
strong winds (except near the windward foot of an iso- 
lated dune). This is what one might have expected in the 
first place, and the reader may wonder at this point how 
the saturation transients and their characteristic length 
scale is can have the claimed importance. How can is 
affect the shape of a dune that is typically about two or- 
ders of magnitude larger? This apparent puzzle is now 
easily resolved by going back to Figs. ||, and Eq.(|25|) 
and by observing that the symmetry breaking shift Sxr 
of the location of the maximum of the shear stress with 
respect to that of the maximum of the height profile (or 
envelope) that is responsible for the finite windward slope 
and growth of dunes, is also of the order of some per cent 
of the total dune length. In summary, the longitudinal 
profile of dunes and heaps is determined by the competi- 
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tion of two quantitatively small but qualitatively crucial 
effects, one related to turbulent wind flow and the other 
to sediment transport. This may be the reason, why its 
explanation proved elusive for a long time. 

To get a qualitative idea of the consequences of the 
introduction of the generalized non-local flux law in 
Eq.(^0|) as replacement for Eq.(||), we want again to go 
back to our discussion of the surface migration velocity 
of the cosine shaped heap in Fig. |8[ Let us for the 
moment adopt a crude approximation and replace the 
expression for given in Eq.(pTl) by its simpler limiting 
form Qs oc T^/^ introduced in Eq.(^. We also neglect 
the variation of £s on the dune and replace it by a (fine- 
tuned) constant £s ~ 0.1 L. With an influx (about 0.7 
for the solid and 0.8 Qs for the dashed line in the lower 
part of Fig. ^) one can thus achieve a fairly constant sur- 
face velocity over the whole length of a cosine shaped 
heap. Again the constancy is slightly better for rt < 2 
than for n = 2. It is further improved by reducing the 
slope of the heap well below the optimum windward slope 
of the dune obtained for — 0, as seen from a compari- 
son of the solid line and the dashed line. We also note in 
passing that the influx needed to maintain the shape is 
increasing with decreasing slope. Together, the plots in 
Fig. 1^ confirms our claim that even an ^ L may visibly 
affect the overall shape of aeolian dunes. Although, for 
the example shown in the lower part of Fig. H, the satura- 
tion length is only about 1/30 of the heap length, the slip 
face instability is evidently completely washed out. Al- 
together, this strongly suggests the existence of transla- 
tion invariant cosine shaped heap solutions for the model. 
The ultimate proof will be provided by the numerical re- 
sults presented in Section M, where the full form of Qs 
and is according to Eq.(|3|) will be used, but the present 
crude approximation already illustrates the main point, 
and also demonstrates that the behavior is a generic con- 
sequence of Eq. (^ and insensitive to the detailed form 
of the parameters is (x) and Qs (x) that may phenomeno- 
logically be somewhat different from the model prediction 
without affecting our general conclusions. 

An immediate consequence of the foregoing discussion 
is the existence of a minimum dune size. For small 
enough dunes, the slip face instability is washed out by 
the saturation transients. One may also arrive at this 
conclusion from an analysis of heaps. To this end we 
observe that the value of the saturation length ig is a 
property of the wind velocity and the saltation kinetics 
and depends on the topography only indirectly through 
the variable shear stress r. Moreover, it is apparent from 
Fig. ^ that this dependence becomes weak for strong 
winds. On the other hand, the symmetry breaking shift 
Sxr is proportional to the absolute size of the heap (or en- 
velope) and not directly dependent on the wind velocity. 
For the special profile /l , this was verified analytically in 
Eq. (p5|) . As we have seen, a smooth heap can only be a 
translation invariant solution of the model if the lag (of 
order is) of q{x) with respect to Qsix) and the shift Sxr 
are fine-tuned to guarantee a vanishing erosion rate at 



the top of the heap. From this we expect heaps to obey 



dxr (X L M const. 



(33) 



to a first approximation. This condition can only be ful- 
filled if the aspect ratio e of heaps grows proportional 
to their size (i.e. roughly e cx H). Hence, in contrast 
to large dunes with slip face, for which we have argued 
that they are asymptotically scale invariant (e const.), 
heaps must have a strongly size dependent aspect ratio. 
As a consequence, translation invariant heap solutions 
obviously cannot exist beyond a certain critical size. A 
slip face will develop when the shear stress on the lee 
side of the heap drops below the threshold value Tf , or at 
the latest, when the lee slope exceeds the critical slope 
for flow separation. This will be further analyzed in Sec- 
tion Finally, we note that the steady-state flux of a 
heap can be estimated by the observation that the outflux 
is essentially determined by the strength of the reduction 
tq — Tniin of the shcar stress at the lee end of the heap. 
According to Eq.(^), the latter is (for a given shape) pro- 
portional to the aspect ratio e. For qualitative purposes, 
the outflux may thus be estimated in the saturated 
flux approximation with Eqs.(31 ),(p3),(P) as 



Qs ^ '^min 



Tf CX £c 



(34) 



where we have assumed Tmin/Tt ^ 2 (fulfilled for mod- 
erate wind speeds and/or heaps near the critical heap 
size) to linearize the expression for qs (r) given in Eq. (^Tj) . 
Here, Sc oc tq — Tf is the critical aspect ratio for which 
the shear stress on the lee drops below the threshold and 
the outflux vanishes. Note that the latter increases with 
increasing shear stress whereas the heap length decreases 
according to Eqs.(^, (^. The effects of the two trends 
onto the critical heap mass could therefore partially can- 
cel unless the lee slope exceeds the critical slope for flow 
separation. 



C. Slip face 

We have argued above that for large heaps (L ^ is), 
aeolian sediment transport tends to increase the down- 
wind slope until it reaches the angle of repose of the 
grains. At this point, any further increase of the lee 
slope initiates avalanches that restore a slope slightly be- 
low the static angle of repose and eventually create a slip 
face of a roughly uniform slope of about 32° — 35°. Since 
the physical modeling of this process itself is not a ma- 
jor objective of the present contribution, we can choose 
between different possible implementations for this phe- 
nomenon. In 2d it is possible to represent the slip face as 
boundary condition for the sand transport. It is uniquely 
determined by its fixed uniform slope and mass conser- 
vation. However, with regard to a future generalization 
of the present 2d model to 3d we chose a more realistic 
implementation based on a widely-used avalanche model 
psf . The formulation of this model bears some close sim- 
ilarities with the sand transport model presented in the 
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FIG. 14: Solution of the minimal model. 



preceding paragraph, and thus suggests itself as a natural 
extension of the latter to the slip face. This completes 
the definition of the minimal model that will be solved 
numerically in the next section. 



V. SOLUTION OF THE MODEL 

Apart from the model definition, the preceding sec- 
tions have provided some qualitative insights into the 
main mechanisms responsible for dune formation and mi- 
gration. Now we are prepared to study numerically the 
quantitative predictions of the model. Again, we em- 
phasize that we only can explore some major features of 
the model in the present report, leaving many interesting 
questions and more systematic and quantitative param- 
eter studies for future work. 

For convenience, the solution procedure of the mini- 
mal model is summarized as a flow chart in Fig. |lj. One 
starts with an initial profile h{x, t = 0) (typically fa or 
Jq), checks whether a separation bubble has to be added 
for the calculat ion o f the shear stress, then obtains the 
latter from Eq.(16a) and uses the result as input for the 
iterative solution of the sand transport equation Eq.(30). 



This finally gives the erosion/deposition needed to up- 
date the surface profile. Technically, Eq.(16a) is imple- 
mented as a Fast-Fourier-Transform algorithm, and for 
the integration of Eqs.(pO|) and an upwind discretiza- 
tion scheme is used. Simulation times can be reduced by 
using an adaptive time step. 



A. Steady— state shapes 

The scheme of Fig. |lj can be iterated for different in- 
flux boundary conditions. For all of the numerical cal- 
culations presented below, we chose periodic boundary 
conditions. They are the natural choice for studies of the 
steady-state shapes. To investigate the mass balance un- 
der prescribed influx conditions, on the other hand, one 
has to apply open boundary conditions. 

Fig. |l5| shows steady-state solutions of the model for 
initial profiles /g of different mass. These solutions 
arc obtained for fixed wind conditions with parameters 
A = 3.2 and B = 0.25 appropriate for the central slice of 
a 3c? (symmetric) heap or of a barchan dune. The shear 
velocity = 0.4 m/s lies well above the impact thresh- 
old. (The situation very close to or below the thresh- 
old would need special attention.) As anticipated above. 




FIG. 15: Steady-state heaps (upper plot) and dunes (lower 
plot). The aspect ratio is stretched for better visualization. 



large dunes become asymptotically scale invariant. The 
asymptotic master curve for the windward profile is prac- 
tically indistinguishable from the profile Jq (n < 2), and 
the slope at the brink is indeed positive. Its average wind- 
ward slope is inversely proportional to the value of the pa- 
rameter A given in Eq .(|l3| ) . Due to the additional terms 
in the expression Eq. (|18|) for the shear stress on dunes 
with a finite width, somewhat steeper average windward 
slopes are predicted for barchan dunes than for trans- 
verse dunes under identical influx and wind conditions. 
However, a detailed quantitative comparison is probably 
beyond the scope of the present semi-quantitative im- 
plementation. More important are the remarkable qual- 
itative predictions of the model. In particular, the fact 
that dunes with slip face are only stable above a certain 
(wind dependent) critical size, whereas smooth steady- 
state heaps only exist below a critical size, deserves at- 
tention. We also note that the steady state is not always 
unique. There is a hysteretic regime, where the initial 
conditions can select one of two possible steady-state 
shapes and accordingly the masses for the two sets of 
profiles in Fig. |l^ are not all distinct. The largest heaps 
in the upper plot were obtained from flat initial profiles 
fa, whereas the smallest dunes with slip face in the lower 
plot were obtained from steeper initial profiles fa of the 
same mass. Especially, the dune with a negative slope at 
the brink could only be obtained from steep initial con- 
ditions. Since under natural wind and sand conditions, 
the initial conditions themselves will generally be heaps 
or dunes close to the steady state, one can say that the 
model predicts a critical heap size for slip face formation 
and a critical dune size for slip face destruction. In both 
cases the slip face is finite as a consequence of flow sep- 
aration. The latter also allows a dune to be somewhat 
higher than a heap of the same mass, since its effective 
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grou nds the simple scaling prediction v oc L ^ in Sec- 
We have also given some arguments why this 
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FIG. 16: Steady-state heights H versus the product of the 
height and length of the heaps and dunes. In the hysteretic 
regime, flat and steep initial conditions have to be distin- 
guished. 




H [m] L [m] 

FIG. 17: Migration velocities predicted by the minimal model 
for steady-state dunes of different size at various wind veloc- 
ities. The caption gives the shear velocity u, in m/s. The 
numerical data are compared to the scaling laws v oc H"^ 
(left) and v oc L~^, where L is the length of the envelope of 
the dune and its separation bubble as described in the main 
text. (Note that the migration of real dunes is substantially 
slower due to the small fraction of wind days per year.) 



volume as seen by the average air flow is increased by the 
separation zone. As anticipated, the aspect ratio of the 
dunes is asymptotically constant, whereas it is strongly 
size-dependent for heaps. This effect can be seen more 
quantitatively in the representation of Fig. |l^, where the 
height H of the steady-state heaps and dunes is plot- 
ted versus the product HL of their height and length 
L. Clearly, heaps are better described by i? oc HL as 
predicted in Eq.(^), whereas large dunes approach the 
scaling limit H oc V HL. 



B. Migration velocity 

For the overall migration velocity of steady-state dunes 
with a scale invariant profile, we derived on general 



tion 

prediction should be rather robust against relaxing the 
condition of shape invariance, in contrast to the relation 

V oc H^^ that can only be inferred from it if the scaling 
assumption holds exactly. Here, we check these predic- 
tions for the steady-state solutions numerically. Fig. |l^ 
shows the numerically obtained migration velocity for 
dunes fitted to the scaling relations v oc H~^ (left) and 

V oc L~^ (right). As we have mentioned above, one has to 
take for L the characteristic length of the envelope rather 
than that of the dune alone. For simplicity, we estimate 
L by adding 6H to the horizontal length from the wind- 
ward foot of the dune to its crest, thus neglecting the 
weak slope dependence of the size of the separation bub- 
ble. Obviously, the L~^— scaling is superior for moderate 
winds and small dunes where v oc H^^ systematically 
fails to describe the data. This is also supported by field 
data 1 22 . Both fits become identical in the scaling limit 



L :$> is- Due to the decrease of £s with the wind speed, 
the latter is reached for smaller dunes at stronger winds. 



C. Stability 

We have already pointed out that the choice of dif- 
ferent boundary conditions for the flux allow a separate 
discussion of shape and mass stability. This is of prac- 
tical importance, since (in 2d) all steady-state shapes 
are unstable with respect to mass changes. If the influx 
of a steady-state solution deviates slightly from its cor- 
responding steady-state flux, this solution will start to 
either shrink until it has flattened out or grow without 
bound. Though the latter effect could (but need not) be a 
peculiarity of the vanishing outflux for 2d dunes with slip 
face, at least the former generalizes to 3d heaps. Despite 
the fact that the steady-state shapes are (locally) sta- 
ble attractors for the shape evolution under periodic flux 
conditions, mass stability can in general not be achieved 
under open boundary conditions. The situation is clar- 
ified in Fig. 18. It depicts the steady-state sand flux q 
over the bedrock as a function of aspect ratio. The nu- 



merical results nicely confirm our theoretical expectation 
from Eq. (^4|) . For all dunes with slip face the flux van- 
ishes identically in 2d, whereas in general it grows with 
decreasing size for smooth heaps. For open boundary 
conditions, the line in Fig. |l^ can be interpreted as an 
unstable phase boundary (with hysteresis) between in- 
finitely growing and shrinking solutions. For example, 
a heap with influx slightly below the steady-state, will 
shrink a bit. To remain close to the steady-state shape, it 
will therefore mainly reduce its height, whereas its length 
will stay almost constant. Due to the reduced aspect 
ratio e, f decreases in magnitude and the shear stress 
depression at the lee boundary is less pronounced. As 
a consequence there is less deposition on the downwind 
slope and the outflux is higher, so that the heap shrinks 
even more etc. A completely analogous reasoning applies 
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FIG. 18: Steady-state outflux under periodic boundary con- 
ditions. In the hysteretic regime, steep and flat initial condi- 
tions have to be distinguished, as in Fig. ^ The figure may 
also be read as a phase diagram for the situation with open 
boundary conditions. In this case the steady-state solutions 
— though attractors for the shape — are unstable with re- 
spect to mass fluctuations. 



to the opposite case of higher influx. The corresponding 
shape attractors are the scale invariant asymptotic dune 
shape and the flat surface, respectively. 

The above discussion explains, why isolated smooth ae- 
olian sand heaps are rarely observed as distinct features 
of desert topographies. Under approximately stationary 
wind and influx conditions they only exist as transient 
states that either vanish or develop into dunes with slip 
face. Under variable wind and influx conditions, the situ- 
ation is less clear and deserves a detailed study of its own. 
For example, the model predicts that during a period of 
strong wind all dunes are driven towards the asymptotic 
shape. After a subsequent period of weak winds, finite 
size effects become more pronounced, and small dunes 
may develop longitudinal profiles like those in the hys- 
teretic regime or even loose their slip face. Again, the 
case To « Tt of a shear stress close to or below the thresh- 
old shear stress needs special attention. 

The prevailing wind conditions as well as recent 
changes in the wind velocity are thus encoded in a compli- 
cated but comprehensible way in the shapes of the dunes 
in a dune field. This is a promising direction for further 
studies. One may hope that by systematic studies along 
these lines one will in the future be able to infer flow con- 
ditions in remote or uncomfortable places (e.g. on the sea 
bottom or on other planets) by analyzing dune shapes. 



D. Relaxation dynamics 

As a first step towards an understanding of the effects 
of variable wind speeds (for constant wind direction), 
this section is devoted to an exploratory study of the 
transient shape evolution. We restrict ourselves to peri- 




b) 

FIG. 19; Growth histories for two Gaussian heaps of different 
mass and initial aspect ratio for periodic flux boundary condi- 
tions. Note that the distances to reach the steady-state shape 
are different. (The aspect ratio of both plots was rescaled by 
a common factor for better visualization.) 



odic boundary conditions leaving the richer phase space 
of open boundary conditions for future studies. Fig. 
shows two extreme scenarios. A flat initial condition with 
a mass greater than the critical mass for slip face forma- 
tion (a), and a steep initial condition with a mass below 
the critical mass for slip face destruction (b). The steep 
initial condition gives rise to the temporary formation 
of a slip face that is finally washed out by the satura- 
tion transients, whereas the flat heap steepens until the 
shear stress on the lee falls to the threshold value. This 
causes complete deposition on the lee side of the heap. 
Whether this happens before or with the onset of flow 
separation depends on wind and influx conditions. As 
it also depends on the precise numerical values of some 
of the phenomenological parameters of the model, a de- 
tailed parametric study is again beyond the scope of the 
present contribution. 

Although the times to reach the steady state are ap- 
parently somewhat longer for the flat initial condition, 
it is evident from Fig. ^ that the relaxation dynamics 
is in general relatively fast even if the initial condition 
is far from the steady state shape. Large dunes under 
low influx conditions as they prevail e.g. in flelds of iso- 
lated dunes should therefore be well described by an adi- 
abatic approximation assuming that (except after drastic 
changes in the wind and sand conditions as they occur 
during sand storms) the dune is practically in a steady 
state. Apart from the low influx, this also relies on the 
fact that virtually no sand is lost over the slip face. For 
a large isolated 3d barchan dune this implies that most 
of the sediment transported over the dune is actually 
trapped in a tread-milling flux, and only a small por- 
tion of the total flux is contributed by and contributes 
to the external flux. Hence, under steady wind condi- 
tions these dunes are in a quasi steady state and thus 
very close to their true steady-state shape. Investigation 
of the steady-state properties is therefore the starting 
point also for the study of their time evolution. More- 



17 



0.2 

A ° 

"-0.2 

H 

_o.4 
-0.6 



0.9 
0.88 
"0.86 
0.84 
0.82 



SXr 

SXa 



0.1 



0.2 



0.3 



0.4 



0.1 



0.2 

time 



0.3 



0.4 



FIG. 20: The figure shows the transient evolution of various 
interesting length scales for a heap. Lower part: Height of the 
heap. Upper part: Distance of the locations of the maximum 
of the shear stress and of the maximum of the sand flux from 
the position of the top of the heap. In the steady state, the 
erosion/deposition vanishes at the crest. 



over, this suggests that a comparison of our steady-state 
shapes to shapes obtained in field measurements is jus- 
tified. In fact, the calculated shapes agree nicely with 
recent measurements for barchan dunes The situa- 
tion is less clear for small heaps, vifhere mass losses can 
be of the order of the total flux and may thus lead to 
significant differences between the steady-state and the 
transient shapes under vanishing influx. 

In the remainder of this section, we want to investi- 
gate more closely the mechanism that drives the shape 
relaxation. As we have pointed out, the positions of the 
maximum of the sand flux and of the maximum of the 
profile must coincide in the steady state to make the ero- 
sion/deposition vanish at the crest. We have shown that 
for small heaps, this can be achieved by a fine-tuning of 
Sxr to about ig. In contrast, for large dunes and strong 
winds, 6xr ^ is, and the steady-state condition can 
only be met with a singularity at the crest. This impor- 
tant difference is exemplified by Figs. EG and 21, Both 



figures show the evolution of the height and the displace- 
ments 6xr and Sxq of the locations of the maximum of 
the shear stress and of the maximum of the sand flux 
from the location of the top of the sand profile, respec- 
tively. The distance between both displacements is the 
lag of the flux with respect to the shear stress due to the 
saturation transients, and is therefore closely related to 
the saturation length is for smooth surface profiles. It 
guarantees the proper vanishing of the erosion rate q' at 
the top of a steady-state heap where 6xt is finite, but 
vanishes for large steady-state dunes, where the slip face 
ends in a sharp brink singularity at which the grains fall 
into the wake and are quickly deposited. 



•7 - 



•-a 




0.1 0.2 0.3 0.4 0.5 0.6 




FIG. 21: The figure shows the transient evolution of various 
interesting length scales for a dune that develops out of a 
smooth heap as in Fig. |l^. Lower part: Height of the dune. 
Upper part: Distance of the locations of the maximum of the 
shear stress and of the maximum sand flux from the position 
of the top of the crest. The lag between shear stress and sand 
flux vanishes, when the slip face reaches the crest. 



VI. SUMMARY AND OUTVIEW 

In summary, we have shown that a simple minimal 
model for the wind-driven sediment transport over a 
sand dune is capable of explaining several important fea- 
tures of desert topographies. Among them are the mi- 
gration velocities of heaps and dunes, their shape along 
the wind direction and the existence of a minimal dune 
size and a maximum heap size. 

As we have emphasized throughout this contribution 
and demonstrated by the numerical solutions in the pre- 
ceding section, the symmetry breaking part of the shear 
stress exerted by turbulent air flow on an obstacle, and 
local deviations of the sediment flux from its equilibrium 
transport capacity ("saturation transients"), are the es- 
sential ingredients in the modeling of aeolian sand dunes. 
It is exactly the balance of these two relatively small ef- 
fects that is responsible for the relaxation of arbitrary ini- 
tial conditions into a characteristic dune or heap shape. 
Their neglect was responsible for the fail ure o f the naive 
zeroth order model discussed in 



Section II B , In hind- 
sight we can say that it is not so much the quantitative 
errors but the omission of this qualitatively important 
mechanism, what makes the zeroth order model an insuf- 
ficient description. In contrast, taking this balance prop- 
erly into account, makes the minimal model structurally 
stable against the neglect of less significant quantitative 
details of the same order of magnitude. 

This direction was recently pursued further in an effort 
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FIG. 22: Qualitative shape diagram that could be useful in 
the analysis and comparison of field studies. The migration 
velocity is constant along rising lines, whereas falling lines 
indicate invariant dune shape. 



to calculate analytically certain steady state shapes of 
dunes and heaps by "linearizing" the minimal model . 
One may as well wish to proceed also into the opposite di- 
rection. After the basic mechanism is understood, more 
elaborate dune models can be constructed by putting 
some of the neglected details back into the description. 
Detailed parametric studies of such a refined model for 
a certain dune type and comparison to field data would 
be very useful to test some of the less generic predic- 
tions of the underlying sand transport model [ic| ], such 
as the shear stress dependence of the saturation length ^s 
(Fig. [T^ ). This is important, since, as we have shown, the 
variable parameter ig sets the characteristic length scale 
with respect to which dunes and heaps can be said to be 
large or small. Phenomcnological knowledge about is is 
still very limited. More detailed studies could further be 
helpful to map out quantitative shape diagrams, of the 
type sketched qualitatively in Fig. |2^. These diagrams 
could be useful not only for the validation of the model, 
but also for the comparison of field data from different 
places with different prevailing wind and sand conditions. 



The migration velocity is constant alon g th e rising lines 
in Fig. |2^, which were obtained from Eq.(|ll|) using qK, 
together with Eq. ( |3l| ) . They allow for example a compar- 
ison of the migration velocities of dunes of different sizes 
that are exposed to (on average) identical winds. Or one 
may infer the average wind speed from measurements of 
sizes and migration velocities in a dune field. The falling 
lines in Fig. ^ are lines of constant shape, assuming that 
the latter is determined by Is/L. They may thus be used 
for correlating wind speeds with dune shapes. In gen- 
eral (in particular for the full 3d problem), such shape 
diagrams will be more complex since the influx is an ad- 
ditional important variable that we have neglected here, 
as it vanishes for 2d dunes in the steady state. 

Moreover, as we pointed out, there are still many con- 
sequences of the present model that await a systematic 
investigation. And a major future task is finally the gen- 
eralization of the present discussion to the 3d case. A 
promising route could be the construction of an effec- 
tively sliced model that allows one to use the proposed 
model for the separation bubble and to keep the time- 
limiting calculation (the integration of the flux equation) 
effectively one-dimensional. The smaller transverse cur- 
rents could be inferred from the sliced solution. A gen- 
eralization of the flux equation to the 2d surface of a Zd 
dune is also feasible A more ambitious task will 

eventually be the simulation of dune fields. Whereas the 
existence of a minimum dune size could be obtained by 
an analysis of the shape stability alone, the question of a 
possible existence of a characteristic or maximum dune 
size in a dune field, depends on the mass balance of a 
dune in the complicated environment provided b y t he 
other dunes, and is much more difficult to answer [ [40| . 
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